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ABSTRACT 

Spitzer-IRS maps of H 2 pure rotational lines from S(0) to S(7) in three out- 
flows from Class sources, namely L1448, BHR71, and NGC2071, are presented. 
These lines are used, in conjunction with available rovibrational, near- infrared 
H2 lines, to probe the physical conditions of the warm gas between hundreds 
and thousands of Kelvin. We have constructed maps of the molecular hydro- 
gen column density, ortho-to-para ratio and volume density, together with the 
index /3 of the power law describing the distribution of gas temperature. In all 
three outflows, the present ortho-to-para ratio significantly deviates from the 
high temperature equilibrium of 3, being on average between 2.0 and 2.3. These 
low values, that in general reflect the young age of these flows, are found also 
in regions of relatively high temperature (around 1000 K), likely indicating that 
shocks are occurring in a time shorter than that needed for a complete para to 
ortho conversion. Density maps indicate upper limits close to LTE conditions, 
i.e. between 10 6 -10 7 cm -3 ; moreover we demonstrate, on the basis of the detec- 
tions of HD emission spots (R(3) and R(4) lines) that a density stratification does 
exist, with the low density components (10 4 -10 5 cm -3 ) associated to the coldest 
gas. The /3 index is found in all flows to be above 3.8: this value is consistent 
with predictions of multiple C-type bow shocks with a range of velocities, some 
of which insufficient to achieve the temperature at which H 2 partially dissociates. 
The contribution of H2 to the total cooling is quantitatively similar to that of 
other abundant molecules emitting in the far-infrared, such as water and CO; 
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moreover, the luminosity radiated away is comparable with the estimated kinetic 
energy of the swept-out outflow: this supports the scenario in which the shock 
from which H 2 is emitted is also capable of accelerating the molecular outflow 
driven at the shock working surface. 

Subject headings: Stars: formation -ISM: individual objects: L1448, BHR71, 
NGC2071- ISM: jets and outflows - infrared: ism: lines and bands 

1. Introduction 

The physical effect of jets from accreting protostars on the parental cloud is the for- 
mation of shock waves. The accelerated, compressed and heated gas radiates away the 
accumulated thermal energy through the emission of lines from various species, in different 
proportions that depend on the relative abundances and the specific shock conditions. In 
this respect, a fundamental role is played by the molecular hydrogen, whose high abundance 
compensates for the fact that its homonuclear nature leads to quadrupole-allowed rovibra- 
tional transitions with small radiative rate. Rovibrational H 2 lines in the near- infrared, and 
in particular the 1-0 S(l) line at 2.12 fim, have been largely used to discover jets from young 
protostars, and, at present around 1000 individual regions of shock emission have been iden- 
tified (Davis et al. 2010). However, the most efficient cooling channels of molecular hydrogen 
are the pure rotational lines of the ground vibrational level, lying in the mid-infrared at A < 
28 jum (e.g. Kaufman & Neufeld 1996). These lines represent a key tool to penetrate the 
very obscured regions that surround the so-called Class sources, i.e. the very young proto- 
stars in the main accretion phase of their evolution (age ~ 10 4 yr, Andre, Ward-Thompson 
& Barsony 1993); moreover the H 2 mid-infrared lines are easily excited at few hundreds of 
Kelvin, and consequently they represent a valuable tool for tracing the cold gas component 
in the jet, where the initial conditions of the emitting gas are often still manifested in a non 
equilibrium value of the H 2 ortho-to-para ratio (Timmermann 1998, Wilgenbus et al. 2000, 
Maret et al. 2009, Barsony et al. 2010). In addition to the different excitation energies, 
the pure rotational lines differ from the rovibrational ones for being characterized by lower 
critical densities: all these circumstances make the combination of the two sets of lines as 
particularly valuable for deriving a global view of gas excitation in molecular outflows and 
for inferring the total H 2 cooling, in particular for the environments of the Class protostars. 
The H 2 pure rotational lines were first observed extensively with the Infrared Space Obser- 
vatory (ISO), by means of the ISOCAM camera and the SWS spectrograph. These studies 
have demonstrated the usefulness of the H 2 rotational lines for deriving fundamental pa- 
rameters of molecular shocks, such as temperature, column density, and ortho-to-para ratio, 
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although averaged on large areas due to the poor spatial resolution of ISO ( Noriega- Crespo, 
Garnavich & Molinari 1998, Neufeld, Melnick & Harwit 1998, Smith, Eisloffel & Davis 1998, 
Nisini et al. 2000, Molinari et al. 2000, Rosenthal, Bertoldi k Drapatz 2000, Benedettini 
et al. 2000, Lefloch et al. 2003). More recently, the IRAC camera (Fazio et al. 2004) 
and the IRS spectrograph (Houck et al. 2004) onboard Spitzer (Werner et al. 2004), have 
allowed us to expand upon the handful of outflows observed with ISO (e.g. Smith et al. 
2006, Cyganowski, et al. 2008, Zhang & Wang 2009, Maret et al. 2009, Ybarra & Lada 
2009, Ybarra et al. 2010, Dionatos et al. 2010, Barsony et al. 2010, Takami et al. 2010); 
more important, however, is that for the first time sensitivity and spatial resolution have 
been achieved close to those reached from ground-based observations in the near-infrared, 
hence allowing H2 observations in different spectral bands to be combined effectively. 
We have started this kind of multiband analysis in a group of outflows from Class sources 
observed with Spitzer-IRS between 5.5 and 37/xm. The data have been presented by Neufeld 
et al. 2009 (hereafter Paper I), who gave an overview of the spectra along with the global 
physical properties averaged over the entire emitting region. The H 2 luminosities were fitted 
by a NLTE model in which an admixture of gas temperatures is present within the spatial 
resolution element. A power-law distribution of the gas temperature is assumed, with the 
column of gas at T and T + dT proportional to T^^dT; (5 typically ranges between 2.3 and 
3.3. Other main parameters of the model are the H2 ortho-to-para ratio, that is found in all 
the flows apart from L1448 lower than the equilibrium value, and the particle density, which 
is typically lower than 10 4 cm -3 . In a second paper, Nisini et al. 2010 (hereafter Paper II) 
have analyzed in deeper detail the outflow of L1157, deriving maps of the main parameters 
(H 2 column density, temperature, ortho-to-para ratio) that govern the shock physics at ~ 
2" resolution, and giving detailed estimates of the molecular cooling. In particular, gradi- 
ents of the H 2 column density are probed along the flow, with the highest values found in 
the blueshifted lobe, where N(H 2 ) ~ 3 1 20 cm~ 2 ; the temperature spans from ~ 250 to ~ 
1500 K, with the largest range of variations towards the intensity peaks close to the exciting 
source; the ortho-to-para ratio is significantly lower than the equilibrium value, especially in 
plateau regions between consecutive intensity peaks. 

The present paper deals with three other outflows of our sample, namely L1448, BHR71, 
and NGC2071, while a dedicated paper will be devoted to the analysis of VLA1623 (Liseau 
et al., in preparation). Together with L1157, they constitute a homogeneous sample of out- 
flows from very young protostars. Each object, however, has its own peculiarity: L1448 is 
a prototype of a source driving a molecular jet (Eisloffel 2000); BHR71 is one of the few 
outflows to be almost fully molecular (Giannini et al. 2004); NGC2071 is an example of 
a parsec-scale outflow of high mass (Stojimirovic, Snell & Narayanan 2008), and L1157 is 
maybe the most active outflow from the chemical point of view (Bachiller & Perez Gutierrez 
1997). 
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The detailed description of morphology and physical characteristics of the three outflows 
discussed here have been given in Paper I and Melnick et al. 2008; here, we limit to briefly 
summarize some information useful for the comprehension of the present paper. The L1448 
outflow (length ~ 0.3 pc), originates from the Class source L1448-mm (L1448C, Bachiller, 
Andre & Cabrit 1991), subsequently identified in Spitzer images as a binary system. At the 
tip of the northern lobe another triple system is located (Looney, Mundy & Welch 2000), 
with a Class source, L1448-IRS3 (L1448-N) driving a compact molecular outflow. Also the 
BHR71 outflow originates from an embedded binary protostellar system (IRS1 and IRS2, 
Bourke 2001), each driving a molecular outflow. IRS1, in particular, has been classified as a 
Class protostar (Bourke et al. 1997); the associated outflow extends in the north-south di- 
rection. The source IRS2 gives rise to a more compact outflow inclined of ~ -36°. Optical [SII] 
emission revealed HH associations (HH320 and HH321) along both the outflows (Corporon 
& Reipurth 1997): these should be at very low degree of excitation, as demonstrated by the 
lack of emission both in Ha and in near- infrared ionic lines (Giannini et al. 2004). NGC2071 
is a parsec scale bipolar outflow originating from a cluster of infrared sources. Among them, 
the source IRS1 has been initially suggested to be the outflow powering source, because it 
dominates the luminosity in the region (Persson et al. 1981); subsequently, however, images 
in H2 emission (Aspin, Sandell & Walther 1992) have demonstrated outflow activity around 
the source IRS3, that is also the major contributor to the luminosity at longer wavelengths 
(Walther, Geballe & Robson 1991), likely due to its young age. 

The main aim of this work is to construct detailed maps of the physical parameters along 
each outflow. Our method is quite similar to that described in Paper II: this will allow us to 
directly compare our results with those obtained in LI 157. Such maps are of fundamental 
importance for elucidating, albeit in a small sample, the typical physical conditions of molec- 
ular outflows, in terms of temperature, density, and mass of the shocked gas. Our results 
also put observational constraints on models for partially (or non)-dissociative shocks. We 
also note that the objects presented here are going to be observed in the far- infrared with the 
Herschel telescope as part of both Guaranteed and Open Time programs. In particular, they 
belong to the sample of the WISH program, devoted to the understanding of the chemistry 
of shocks in protostellar outflows: the maps and the physical parameters given in the present 
paper will be of fundamental importance in deriving the water abundance along the flow, 
for example. 

The paper is organised as follows: in §2] we briefly describe the observations and the results 
obtained; §3] is dedicated to the analysis of the H2 emission, in which we have adopted both 
an LTE and non-LTE approximation for the line analysis, while in §3]we analyze a few spots 
where HD lines were detected. §21 and §51 are dedicated to a brief comparison with studies 
of Paper I /II and a discussion on the outflows' energetics; finally, in §7] we summarize our 
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conclusions. 



2. Observations and results 

All the three outflows have been observed during the Cycle 4 of the Spitzer mission 
(February-March 2008), with the exception of a single V x 1' region mapped previously 
toward each of BHR71 and NGC2071. IRS observations in the full spectral range (5.2-36.5 
/xm) were carried out, with the Long-High (LH), Short-High (SH) (A/AA ~ 600) and Short- 
Low (A/ A A between 64 and 128) modules. The total integration time in each object is 
between 19 and 23 hours. Each outflow map has been obtained by combining several V x 1' 
maps, which were arranged along the outflow axis by stepping the IRS slit by one-half of 
its width in the direction perpendicular to its length. In the SH and LH modules, the slit 
was also stepped parallel to its length by 4/5 (SH) and 1/5 (LH). Further details on the 
observational strategy along with the data reduction procedure are given in Paper I. 
The molecular emission is largely dominated by H 2 emission that has been detected in all of 
three outflows. All the pure rotational lines of the ground vibrational level accessible with 
Spitzer, i.e. from S(0) to S(7), have been detected. In addition, a few water lines (whose 
analysis will be the subject of a future paper), and the HD R(3) and R(4) rotational lines 
have been detected in a few spots along the L1448 and NGC2071 outflows. In Table 1 we 
report the H 2 and HD brightness in selected areas (of radius 5") that represent peaks of 
emission or mid-infrared counterparts of peaks at other wavelengths. In the same table, we 
also give the brightness of the 1-0 S(l) line at 2.12 /xm computed for the same area. The 
1-0 S(l) image of BHR71 was retrieved from the ESO archivd^, while the other two were 
taken from Davis & Smith (1995) and Eisloffel (2000). The images were first rebinned at the 
IRS spatial sampling (of ~ 2"/ px) and then flux-calibrated by using bright field stars whose 
photometry was retrieved from the 2MASS catalog. 

Figures [U [2] and [3] show the images of the 0-0 S(l) and 0-0 S(5) lines, overlaid with contours 
of the 2.12 /xm line intensity in steps of 20% from the peak. The overall morphology of the 
0-0 S(l) and 0-0 S(5) lines is similar; this latter, in particular, correlates well with the peaks 
of the 1-0 S(l) line (right panels), given the quite similar excitation energy, that also makes 
them probes of the same gas component at thousands of Kelvin. 

Figure [1] shows the line maps of L1448, where we have labelled the peaks of the H 2 emission 
along with the positions of the two HD peaks that are revealed by the R(3) and R(4) lines. We 
have identified two peaks (N and S), in the northern and in the southern lobe, respectively; 
peak N, in particular, is few arcseconds south-east of the source L1448-IRS3. Following the 
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Fig. 1. — L1448 maps of the 0-0 S(l) and 0-0 S(5) lines, this latter with overlaid contours of 
the 1-0 S(l) line in steps of 20% of the peak, which is 5.0 10~ 5 erg s _1 cm -2 sr _1 . Vertical 
bars indicate the 0-0 lines brightness in the same units. Regions labelled HD1 and HD2 
represent peaks of HD emission, where the R(3) and R(4) are detected, while regions N and 
S are peaks of H2 emission in the northern and southern lobe, respectively. Following the 
nomenclature of Dionatos et si. (2009), we have also marked the 'OF2' region, where other 
near-infrared H 2 lines have been detected. 
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Fig. 2. — As for in FigureHJfor BHR71. The peak of the 1-0 S(l) line is 1.8 10~ 4 erg s _1 cm -2 
sr _1 . The region south IRS1, where the 0-0S(5) is not detected (in contrast with the 1-0S(1) 
line) was contaminated by strong stripes and artifacts, that have been removed within the 
data reduction procedure. The positions of two Herbig-Haro objects (HH 320/HH 321) and 
that of a peak of SiO emission (Gusdorf et al. 2011), are indicated. 
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Fig. 3.— As for in Figure Q] for NGC2071. The peak of the 1-0 S(l) line is 9.0 10~ 4 erg 
s _1 cm" 2 sr -1 . The central region close to IRS 1 results saturated in the IRS map and the 
very southern bow-shock was not mapped. Regions labelled HD1 and HD2 represent peaks 
of HD emission, where the R(3) and R(4) are detected. 
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nomenclature of Dionatos et al. (2009), we have also labelled the 'OF2' position, for which 
near-infrared spectra have been obtained, and that we will discuss in greater detail in §3.2.21 

In Figure [H we have labelled the mid-infrared counterparts of the two Herbig-Haro 
objects in BHR71 (HH320 and HH321), firstly seen by Bourke (2001) at 2.12/im and then 
observed spectroscopically in the near-infrared by Giannini et al. (2004). In the same image, 
we have also identified the mid-IR counterpart of a SiO peak located in the northern lobe of 
the outflow excited by the source IRS1, which has been recently investigated by Gusdorf et 
al. (2011). In agreement with that paper, we label this peak as 'SiO knot'. 
Finally, Figure [3] represents the portion of the outflow of NGC2071 mapped with IRS. The 
overall morphology of the S(5) line follows that of the 1-0 S(l) line, although the IRS maps do 
not cover the very southern region and are saturated around the position of the source IRS1. 
We detected two spots of HD emission, one in each outflow lobe, that are also correlated 
with the H 2 peaks. 

3. Analysis of the H 2 emission 

3.1. Maps of averaged parameters 

In this section, we describe the analysis of the H2 lines, used to derive the maps of the 
main physical parameters characterizing the molecular gas. As a first step, we will analyse 
the pure rotational lines only: given their relatively low excitation energy, most of them 
are observed both along the outflow axis and in the external parts where the gas interacts 
with the environment: consequently, we will be able to probe portions of the outflows where 
gradients in the parameters are expected to occur. As a first step, we assume LTE conditions 
for the 0-0 lines. This assumption is supported by the low values of the critical densities, 
that never exceed 5 10 6 cm -3 for T between 500 and 2000 K (Le Bourlot, Pineau des Forets 
& Flower 1999). For the analysis, we consider only collisions with H 2 molecules, since the 
critical densities do not substantially change if we include also collisions with H and He. 
Possible departures from LTE will be considered further on (see $3j2). Our method follows 
that of Paper I and Paper II: for each pixel we have constructed the rotational diagra nfl from 
the S(0) to S(7) lines, plotting the column density in each rotational state as a function of 
the excitation energy. Such a diagram resembles a straight line if all the emitting gas is in 
LTE at the same temperature (e.g. Gredel 1994), while it follows a curved distribution in 



2 This is the plot ln(N„j/gjg s ) vs. E„j, N^j being the column density of the level with vibrational state 
v and rotational state J, s the spin quantum number, g the degeneracy, and E^j the excitation energy of 
the upper level (in K). 
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the more general case where gas at different temperatures are within a spatial resolution 
element and/or are intercepted along the line of sight. Since we are observing post-shocked 
regions, where the temperature ranges from thousands to hundreds of Kelvin, we have fitted 
the observations by assuming a power law temperature distribution of the form dN oc T^ 13 
dT. Shock models predict (3 typically > 3.8 in bow C-type shocks (Neufeld & Yuan 2008), 
while lower values are expected in J-type shocks, where the shock conditions are strong 
enough to excite lines at high excitation energy (up to ~ 30 000 K, e.g. Flower et al. 2003). 
Consequently, we have taken (3 as a free parameter. It is important to underline that also 
a change in the H 2 density influences the curvature shape in the rotational diagram due to 
the different critical densities of the lines. This implies that a certain degree of degeneracy 
exists between f3 and n(H2), and that, for the LTE conditions we are assuming, the (3 values 
we fit have to be considered as upper limits. 

To derive the H 2 column density, we have integrated the power law between a minimum 
(T m j n ) and maximum (T max ) value of temperature. Since we expect that temperatures 
at thousands of Kelvin contribute to the column density of the pure rotational lines only 
marginally, we have fixed T ma x at 4000 K (and we have verified a posteriori that N(H 2 ) 
varies of only few percent if T ma x = 5000 K). Conversely, since we expect that the bulk 
of the emission traced by the observed lines comes from cold gas, we have chosen not to 
take T m in as constant over the entire flow; rather we adopted for each pixel as T min the 
temperature fitted in the rotational diagram of the S(0) and S(l) lines only, i.e. the coldest 
gas component we are able to trace. T min ranges typically between ~ 150 and 400 K. 

Another free parameter in our fit is the ortho-to-para (OPR) ratio value, which can sig- 
nificantly depart from the equilibrium value of 3 at the low temperatures probed by the 0-0 
lines. Non-equilibrium values are evident in a 'zigzag' behaviour of the rotational diagram 
(see Figure 16 of Paper I), in which the states with even- J lines lie systematically above those 
with odd- J. The departure from the equilibrium value decreases gradually with increasing 
J, such that the highest energy levels (corresponding to S(6) and S(7) lines) have an OPR 
~ 3; we have however considered a single value for all the lines in order to minimize the 
number of free parameters. This assumption thus implies that the OPR we derive is indeed 
an average value. 

Finally, the other parameter affecting the observed line intensities, and consequently 
the column densities, is the visual extinction. Since this affects very marginally the fluxes of 
the 0-0 lines, we have taken Ay constant toward each shocked region. The extinction values 
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have been taken from literature: we have assumed Ay=2 mag for BHR71 (Giannini et al. 
2004), A y =6 mag for L1448 (Nisini et al. 2000), and A v =13 mag for NGC2071 (Melnick et 
al. 2008). In the range from 1 to 13 fim we adopted the extinction law of Rieke Sz Lebofsky 
(1985) and that of Mathis (1998) from 13 to 23 /im, which was extrapolated at 28 /im. 

Summarizing, our fits have three free parameters: index of the dN vs. T power law, 
the OPR, and the total column density N(H 2 ). These three parameters have been fitted for 
each pixel where at least four lines are detected at S/N > 3 by means of a x 2 minimization. 
Typical uncertainties in the parameters are as follows: N(H 2 ): ~ 10%, OPR: ± 0.2, 0: ± 0.2. 
Together with the maps of the above parameters, we also derived two maps of temperature, 
that represent an estimate of the maximum (T warm ) and the minimum (T co id) temperatures 
contributing to the observed emission; these have been obtained as fits of the S(0)-S(l)-S(2), 
and the S(5)-S(6)-S(7) lines, respectively, once corrected for the value of the OPR in the 
considered pixel. 

In Table [2] we give the value of each parameter averaged over the entire map, with, in 
parentheses, the minimum and maximum value attained at particular positions along the 
flow. At this stage we do not give the fitted values of 0, which will be derived again in NLTE 
approximation, i.e. once the degeneracy with n(H2) has been removed, see §3.2.11 
From an inspection of Table 2, it is clear that all the parameters vary significantly, and 
that similarities and differences are apparent among the objects. These are discussed below 
with reference to the individual maps (Figures H]- E]). In the top left panel we present the 
H 2 total column density: a good correlation between N(H 2 ) and the intensity peaks in the 
S(0) line (whose contours are overlaid) is clearly evident, implying that the large majority 
of the molecular gas is at low temperature. Variations up a factor ~ 4-7 are evident in all 
the three flows, with N(H 2 ) ranging between 0.4-4 10 20 cm -2 in L1448, 3-20 10 19 cm -2 in 
BHR71, and 5-18 10 20 cm" 2 in NGC2071. In L1448 (Figgj), the maximum N(H 2 ) is attained 
close to the apex of the northern bow shock: this is the region where the gas is strongly 
compressed because of the interaction of the high- velocity outflowing gas with both the high- 
density clump hosting the system of L1448-IRS3, and the associated outflow; more constant 
values occur in the southern, red-shifted lobe. We do not exclude, however, the possibility 
that a significant column density variation can occur at the apex of the bow, which was 
not observed with IRS. A similar behaviour is recognizable in the N(H 2 ) column density 
map of BHR71 (Figj5]), where peaks of more than 10 20 cm -2 are found at the apex of the 
two outflows driven by sources IRS1 and IRS2. These are the regions where stronger shock 
conditions are evident both at optical-infrared wavelengths (HH 320/321) and in the sub-mm 
(SiO knot). A more constant behaviour is found in the N(H 2 ) values fitted in the NGC2071 
outflow (Figj6]). Both lobes present peaks of N(H 2 ) ~ 10 21 cm -2 , particularly in the parts 
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closest to the cluster of young sources at the centre of the two lobes. 

The OPR, shown in the top right panels of Figures HHHl presents significant departures from 
the equilibrium value of 3 in all three outflows. Similar values have been found, e.g., on 
NGC1333 (Maret et al. 2009), on the Wasp- Waist nebula (Barsony et al. 2010), and on 
L1157 (Paper II). In particular, low values of OPR, typically between 2.0 and 2.5 are often 
found along the external parts of the shock, likely because these regions are where the warm 
gas in the jet entrains the cold one of the ambient medium. Conversely, the highest values 
are generally found where gas components up to ~ 2000 K are probed (see below); we note, 
however, that the equilibrium value of 3 is never reached: this indicates that the shock 
activity, presently on-going, heats the gas for timescales shorter than those required for a 
complete para-to-ortho conversion. 

In the bottom right panels of Figs. HUH] we show the maps of T mrm . As anticipated above, we 
derive this quantity from the linear fit through the most excited 0-0 lines, i.e. S(5)-S(6)-S(7). 
These are detected in L1448 and in NGC2071 only close to the jet axis, while in BHR71 the 
high-J emission appears more extended. T warm ranges typically from 1000 to 2000 K, with 
lower and almost constant values found along the jet axis of L1448 and NGC2071; in BHR71 
peaks of temperatures are found at the locations of the HH objects and close to the SiO 
knot. The typical difference between T warm and T co id (see bottom left panels), which gives 
an idea of the temperature gradient, is between 900 and 1200 K: this indicates that these 
'warm' regions are also dynamically active. Of course, we cannot derive a similar conclusion 
on the regions where the high-J lines remain undetected; however we note that these are the 
same where also T co ^ is low (typically ~ 200-250 K): hence it is reasonable that, even if a 
gradient up ~ 1000 K should exist, the maximum temperature would not be high enough to 
excite bright high-J lines. 



3.2. NLTE analysis 

In this section, we will relax the assumption of the LTE, valid for the pure rotational 
lines, by including in our analysis near-IR rovibrational lines coming from levels with critical 
density significantly higher than the v =0, J< 9 levels. In particular, maps of the 2.12 /im lines 
will be used to infer on the H 2 density (SJ3J2.1), while spectroscopic observations performed 
in selected regions will be used to construct rotational diagrams covering a large range of 
excitation energies (§212.2). 
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Table 1. H 2 photometry in selected areas. 
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a Brightncss computed over a 5" FHWM Gaussian aperture around the indicated position. Ab- 
solute uncertainties are between 20% and 40%. 

b Positions respect to L1448-mm (3h 25m 38.9s 30d 44'06"). 

Positions respect to IRS1 (12h 01m 36.7s -65d 08'48".5). 

d Positions respect to IRS1 (5h 47m 4.8s Od 21'42".8). 
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Fig. 4. — Map parameters of L1448. Top panels: total H2 column density (left); ortho-to-para 
ratio (right); bottom panels: 'cold' (left) and 'warm' (right) gas temperature components, 
respectively derived from the fit through S(0)-S(l)-S(2) and S(5)-S(6)-S(7) lines. In the top 
left panel contours of brightness of the S(0) line are shown for comparison, plotted in steps 
of 20% from the peak. Vertical bars indicate the values of the parameters in the units given 
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- 16 - 




100 -100 100 -100 

Offset (arcsec) 



Fig. 6.— As in Figure 1 for NGC2071. 
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3.2.1. NLTE maps 

As anticipated above, we will exploit the high critical density of the 2.12 /im line (ranging 
from 1.8 10 10 and 4.2 10 8 cm" 3 for 1000 K < T < 2000 K and collisions with H 2 only), to fit 
the lines under NLTE conditions and to derive a map of the H 2 particle density. Moreover, 
since we are deriving n(H 2 ) independently from the Spitzer lines, we also expect to remove 
the degeneracy between n(B.2) and (3 in the rotational diagrams. We stress that all these 
parameters (as shock velocities, temperature, and ionization), are expected to change rapidly 
on sub-arcsecond scales, ie. much smaller than our spatial resolution. In this sense, the 
densities and (3 values derived here are therefore mean values. Our NLTE model consists 
of 24 and 26 levels for ortho- and para-H 2 , respectively. We have considered collisions with 
H 2 and atomic hydrogen, whose de-excitation coefficients are taken from Le Bourlot, Pineau 
des Forets & Flower (1999). We have taken as input parameters the values of N(H 2 ), T min , 
and OPR fitted for each pixel in LTE approximation. Indeed, the inclusion in the fit of 
the 1-0 S(l) line does not affect either the first two parameters, which depend mainly on 
the column density of the 0-0 lines at low- J, or the value of the OPR, which was obtained 
as an average over all the 0-0 lines, some of which very close in excitation energy to the 
level emitting the 2.12 /zm line. Hence, we are left with three free parameters: n(H2), (3 and 
the fraction of hydrogen nuclei in atomic form, / = n(H)/[2n(H 2 )+n(H)]. In particular, as 
shown in Paper II in the outflow of L1157, collisions with atomic hydrogen can indeed play 
an important role in partially dissociative shocks, given their high efficiency in exciting H 2 . 
To evaluate the effect of atomic collisions, we have considered two cases, fully molecular gas 
(/=0) and part of hydrogen in atomic form (/=0.5//=0.9); the results for f=0 are shown 
in Figures EE 

In L1448 (Figure [7]), we generally find higher densities in the nothern lobe, with a density 
peak close to the outflow tip and to source IRS3; in this part of the outflow, the gas appears 
close to thermal equilibrium, with values of n(B.2) close to 10 7 cm -3 . By contrast, in the 
southern lobe n(Hi) is more constant and never exceeds ~ 1-2 10 6 cm~ 3 . 

The effect of increasing the assumed value of / is to decrease the derived H 2 density 
(Paper II). We have verified that a decrease by a factor between 2 and 5 is found all over the 
map if /=0.5 and up to a factor of 10 if /=0.9. Hence, n(B.2) values shown in Figure [7] have 
to be considered as upper limits. 

In the right panel of Figure [7J we show the map of the index 0. We fit values typically 
between 3.5 and 4.5. These values are those expected in bow-shocks where a small fraction 
of H 2 is dissociated: indeed a value of ~ 3.8 is predicted for a single parabolic bow-shock 
where the apex velocity is above the limit for H 2 dissociation; yet higher values are expected 
for a mixture of bow-shocks, including fully non-dissociative shocks. These latter, indeed, 
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are fast enough to excite only the low-lying H2 levels, therefore increasing the slope in the 
rotational diagram of the function dN vs, T (Neufeld & Yuan 2008). We find values of 
that progressively become higher with increasing distance from the central source: this could 
result from an increasing number of unresolved bow-shocks with different apex velocities. 
This result is supported by echelle spectroscopic observations (Davis & Smith 1996) which 
evidence very wide (80-100 km s" 1 ), typically bow-shaped profiles. 

The density map of BHR71 is presented in Figure [SJ A spread of about an order of 
magnitude, from few x 10 6 cm -3 to more than 10 7 cm -3 is found; in general, larger densi- 
ties are found in the southern lobe, even if some enhancements up to 10 7 cm -3 are found 
at the SiO knot and in HH320. As in L1448, we obtain a decrease in n(H2) of up a fac- 
tor of 5/10 if /=0.5/0.9 is taken. The index 0, presented in the right panel, is generally 
lower in the northern lobe, where values exceeding 4.0 are more rarely found; in the southern 
lobe, ranges between 3.5 and 4.5, with the highest values found at the apex shocked region. 

NGC2071 (Figure [9]) is the object where we find the highest densities, close to LTE along 
both the outflow lobes; indeed densities higher than 2 10 6 cm -3 are derived also if /=0.9 is 
assumed. The estimated value of is always larger than 4, with a patchy behaviour that 
argues for the presence of multiple unresolved bow shocks along the line of sight and/or in 
the spatial resolution element. 

3.2.2. Analysis of H2 ro-vibrational lines in selected regions 

In the previous setions we analysed 2-D maps of the shock parameters, basing on the 
Spitzer and 1-0 S(l) line images. Now we will focus on few selected regions, i.e. the 'OF2' 
position in L1448 (Dionatos et al. 2009) and HH320 and HH321 in BHR71, labelled in 
Figured] and Figure[2j In these regions a number of H 2 ro-vibrational lines have been detected 
by means of low-resolution spectroscopy between 1.0-2.5 /mi (Dionatos et al. 2009, Giannini 
et al.2004). We remark that, in addition to further constrain the physical conditions of these 
three particular regions, the inclusion in the rotational diagram of other NIR lines gives us 
the opportunity to check how much the determination of the hydrogen density depends on 
the number of the involved lines and hence to validate (or not) the approach described in 
the previous section, where maps of n(H 2 ) were constructed relying on the 1-0 S(l) line only. 

The inter-calibration factor between Spitzer images and long-slit data has been found 
by scaling the line flux of the 2.12 /im line measured in the slit to the photometry on the 
2.12 /im image computed on the same areas as in the Spitzer images (see Tabled]). The same 
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Table 2. Parameter statistics". 



Object 


N(H 2 ) 6 


M(T>T mm ) b 


OPR 6 


^p6 

cold 


warm 


/3 C 


n(H 2 ) c - d 




(10 20 cm" 2 ) 


(10- 2 M ) 




(K) 


(K) 




(10 6 cm" 3 ) 


L1448 


1.1 (0.4-4) 


1.0 


2.3 (0.5-2.8) 


310 (270-440) 


1250 (1000-1300) 


4.1(3.5-4.5) 


0.3-2 


BHR71 


1.0 (0.3-2) 


0.6 


2.0 (1.7-3.0) 


310 (250-400) 


1500 (1100-1650) 


3.8(3.0-4.5) 


2-7.5 


NGC2071 


7.5 (5-18) 


26.2 


2.2 (1.9-2.9) 


240 (180-300) 


1250 (1100-1400) 


4.2(3.6-4.5) 


6-19 


L1157 


1.4 (0.5-3) 




1.8 (0.4-3.0) 


280 (220-470) 


1230 (1000-1550) 


4.5-5. 5 e 





a Paramcter value averaged over the entire map along with its range of variation (in parenthesis). 
b Parameter derived in LTE approximation. 
c Parameter derived in NLTE approximation. 
d Values referring to f=0.5 and f=0, respectively. 
°Obtained under LTE approximation by N10. 




Offset (arcsec) 

Fig. 7. — Map parameters obtained under NLTE in L1448. Left: map of H 2 particle density; 
no collisions with atomic hydrogen are considered; right: map of the power low index (3 of 
the relationship dN oc dT. 



-20 - 




Fig. 8.— As in Figure for BHR71. 
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Fig. 9.— As in Figure for NGC2071. 
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factor has been applied to the fluxes of other NIR lines, by assuming the line ratios with the 
2.12 /im measured with the spectroscopy. For HH320 and HH321, we have taken an average 
of the fluxes measured over the sub-structures (HH320A/B and HH321A/B) that appear 
resolved at the spatial resolution of the long-slit observations. We have considered lines with 
S/N larger than 5; this implies that lines with upper vibrational quantum number (v up ) up 
to 2 and 3 are considered for L1448 (0F2) and BHR71 (HH 320/321), respectively. 
With respect to the NLTE approach described in the previous section, a further modification 
has to be introduced when lines with higher vibrational number are included in the model. 
Indeed these lines probe gas at temperatures significantly higher than the pure rotational 
lines, at which the OPR is expected to have reached the equilibrium value of 3. Indeed, 
no deviations from this value are recognizable in the rotational diagrams of the NIR lines 
in HH320 and HH321 (Giannini et al. 2004). Consequently, a para-to-ortho conversion 
time as a function of temperature has to be included in the fitting procedure, so as to 
allow lines excited at different temperatures to have different OPRs. Following the approach 
described in Paper I and Paper II, we have considered the variation of the OPR as a function 
of the temperature from a starting value OPRo for a gas that was originally heated at a 
temperature T for a time r. Assuming that the efficiency of the para-to-ortho conversion is 
regulated by the efficiency of the collisions with atomic hydrogen, the dependence of OPR 
with temperature can be expressed as: 



OPR(t) _ OPRp - n (H)kT , OPRlTE /j _ - n {H)kr\ m\ 

1 + OPR{r) 1 + OPR 1 + OPR L te [ ' 

In this expression, n(H) is the number density of atomic hydrogen and OPRlte is the 
ortho-to-para ratio equilibrium value. The parameter k is given by the sum of the rates 
coefficients for para-to-ortho conversion (k po ), estimated as 8xlO _11 exp(-3900/T) cm 3 s -1 , 
and for ortho-to-para conversion k op ~ k po /3 (Schofield 1967). Thus the dependence of the 
OPR on the temperature is implicitly given by the dependence on T of the k coefficient. In 
practice, the inclusion of the NIR lines in our fit introduces as further free parameters both 
the coefficient K— n(H)r and the initial value OPRo- 

In Figure [TOj we show the best fit we obtain for the HH objects in BHR71. First, we note 
that n(H 2 ) differs by about an order of magnitude between HH320 and HH321, with values 
within a factor two with respect to those reported in Figure El From the ratio N(H2)/n(H2) 
we can estimate the H2 cooling length: this is ~ 2 10 13 cm and 5 10 12 cm for HH320 and 
321, respectively. Both these determinations are well below the cooling lengths predicted 
by non-dissociative, C-type shock models, which are never less than 10 15 cm (Neufeld et al. 
2006), while are in better agreement with predictions for slow, dissociative shocks (v s h oc k ~ 
20-25 km s -1 , Wilgenbus et al. 2000). As far as /3 is concerned, and again with reference 
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to Figure El we find a good agreement in HH320 and a value a bit higher in HH321: we 
ascribe this discrepancy to the poorer quality of the fit (with respect to HH 320) . Other two 
parameters (N(H 2 ) and T min ), are fully consistent with those originally found by assuming 
LTE. In regard to the newly determined parameters, we consistently find initial values of 
the OPR of 1.2 and 1.4 in HH320 and in HH321, respectively, i.e. significantly lower than 
the values between 2 and 2.5 of the present OPRs in the same positions (see Figure [8]). The 
parameter K= n(H)r is fitted to the value 1 10 6 yr cm -3 for both the HH objects. Inter- 
estingly, we can remove the degeneracy between n(H) and r, by assuming for this latter 
the age of the shocljfl determined by Gusdorf et al.(2011); the best fit gives an age of 815 
yr and 700 yr for HH320 and HH321, respectively. These imply, for HH320, n(H) = 1.2 
10 3 cm" 3 , and / ~ 4 10" 4 , and, for HH321 n(H) = 1.4 10 3 cm" 3 , and / ~ 5 10~ 5 . Hence, 
we conclude that most of H 2 remains in molecular form. This result agrees both with the 
non-detection of Ha (Corporon & Reipurt 1997) and [Fell] emission (Giannini et al. 2004), 
and also with the results coming from shock models, which are consistent with low-velocity, 
partially dissociative shocks (Gusdorf et al., 2011). Interestingly, significantly higher values 
of / (~ 0.3) are found in Paper II in L1157, where evidence for a large fraction of ionized 
gas has been recognized (Caratti o Garatti et al. 2006) 

In Figure [HI we present the rotational diagram for the OF2 position in L1448. Although 
fewer rovibrational lines are observed here, we successfully constrain all the parameters in 
the fit, finding a good agreement with the values from the LTE model. In this case, we find 
a very low value of OPRo (0.9). For comparison, initial conditions very close to those of 
the cold molecular cloud, where H 2 is expected to be almost all in the para form, have been 
invoked by Maret et al. (2009) to explain the extremely low inital OPR (~ 0.5) measured 
toward NGC1333. Conversely, the higher values found in BHR71 could indicate that inter- 
mittent shock events have occurred, with a temporal interval shorter than the time needed 
to H 2 to cool enough to convert all the ortho states in the para form. 



4. Analysis of the HD emission 

In L1448 and in NGC2071, we detected two HD lines, the 0-0 R(3) transition at 28.50 
/jm and the 0-0 R(4) transition at 23.03 jum, in two distinct positions along the outflows 



3 The shock age is defined as the ratio between the length of the shock wave (roughly comparable with 
the dimension the shocked region) , and the shock velocity. 
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open circles indicate observational and fitted data, respectively; lines coming from different 
vibrational states are indicated with different colours. The parameters of the best fit through 
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as in Figure [TTJ1 
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(see Figures [I] and |3]) . The observation of HD lines in NGC2071 was presented previously 
by Melnick et al. (2008), who used them to infer a fractional abundance of HD of ~ 1-2 
10" 5 . Here, we will use the HD detections to derive the physical conditions of the emitting 
gas independently from the H 2 emission. 

The R(3) and R(4) lines originate from levels with excitation energy of 532 and 883 cm -1 , 
lower than the level v=0, J=4 (1168 cm -1 ) emitting the S(2) line. Hence it is reasonable 
to assume that the bulk of the observed HD emission originates from the same (cold) gas 
component that emits also the H 2 S(0)-S(2) lines. The closeness of the excitation energy of 
the two HD lines makes them not particurlarly suitable for determining the temperature; 
however, unlike the H 2 0-0 lines, they are characterized by higher critical densities (around 
5 10 4 - 10 5 cm~ 3 at 200 K) and hence they can efficiently probe the molecular gas density. 
To this end, we have applied a 8 level NLTE model for HD emission, including collisions 
with H 2 and atomic hydrogen. Collisional and radiative coefficients were taken from Flower 
et al. (2000). As output, we constructed the diagnostic diagram shown in Figure [121 We 
plot the ratio of the HD lines R(3)/R(4) against that of the H 2 lines S(2)/S(0) over a grid of 
temperature (from 200 K to 400 K, continuous lines) and molecular hydrogen density (from 
3 10 4 to 10 6 cm -3 , dotted lines). With a long dashed line we show how the R(3)/R(4) ratio 
changes in the case n(H. 2 )=?> 10 4 cm" 3 and / = 0.5. Superimposed on the diagram are the 
ratios measured in the HD spots of L1448 and NGC2071 (labelled in Figured] and Figure|3]). 

From this diagram we note that: 

- the density is lower than few 10 5 cm -3 in NGC2071 and 10 6 cm -3 in L1448. This is 
in evident contrast with the maps we have derived with the 2.12 /xm line, where n(H 2 ) 
ranges between 10 6 and more than 10 7 cm" 3 . The two results cannot be reconciled 
even if a very large fraction of atomic hydrogen is assumed, since the effect would be to 
further decrease n(H 2 ) (see also §3. 2. H and the caption of Figure [12]). More likely, the 
cold gas component, namely that traced by the low- J H 2 lines, is at density lower than 
the warm gas traced by the rovibrational lines. This implies that, on the one hand, 
a density stratification of the molecular gas does occur, and, on the other hand, that 
the densities derived with the 2.12 fim have to be considered as an average over all the 
temperature components, and as stringent upper limits for the cold gas component. 

- interestingly, the data obtained at HD2 in L1448 indicate a gas temperature lower than 
200 K. This result comes directly from the absence of lines with J > 1 in this position, 
as evidenced in FigurefTB"! We estimated the H 2 column density at HD2 by considering 
the detection of the S(0) and the 3 a upper limit on the S(l) line. We find N(H 2 ) > 5 
10 20 cm -2 , larger than other estimates along the L1448 flow (Figure H]): hence we can 
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explain the detection of HD lines at the HD2 location with a local enhancement of the 
column density of the coldest gas component close to L1448-mm. 

- in the other three spots of HD (HD1 in L1448 and HD1/2 in NGC2071), the gas tem- 
perature exceeds 250 K: according to our NLTE model these conditions are favourable 
for having bright HD lines. Possible local enhancements of the HD abundance can also 
explain the detection of the HD spots: this topic will be treated in dedicated paper 
(Yuan et al. in preparation). 



5. Comparison with results of Paper I and II 

The first aim of this section is the comparison of our fit results with those obtained 
on the same sources in Paper I, which were however obtained by averaging the observed line 
intensities all over the flow. In that work, the H 2 rotational lines have been fitted in a NLTE 
framework, taking as free parameters M(H 2 ), (3, n(H 2 ), OPRo, and n{H)r. With reference 
to Tabled of the present work and Table 3 of Paper I, we first note that marginal differences 
exist among N(H 2 ) and M{H 2 ) derived in the LTE and NLTE treatments. In general we 
find sligthly higher values (within a factor of two), likely arising from having included in our 
fits the S(0) line, which probes the emission of the diffuse molecular gas. In the case of the 
parameters we determine in NLTE, i.e. (3 and n(H 2 ), the two approches give substantally 
different results. The best fits of Paper I give (3 between 2.3 and 3.3 and log(n(H 2 ))~ 3.8. By 
contrast, we get (3 always larger than 3.8, and \og(n(H 2 ))> 6; moreover, though at a lesser 
extent, n(H 2 ) remains larger than the Paper I determinations even in locations where we are 
able to remove the degeneracy with (3. We ascribe this major discrepancy to the sensitivity 
of the 0-0 lines only to the low-density components, a problem that has been overcome in 
the present work thanks to the inclusion of the 2.12 /im line in the analysis. Noticeably, for 
the same reason, in the northern lobe of NGC2071 we probe densities at least an order of 
magnitude higher than those estimated by Melnick et al. (2008) with the Spitzer lines. 

Finally, it is interesting to compare the parameters of the three outflows considered here 
with those obtained in Paper II on L1157. These are summarized in the last row of Table 
[2j The values fitted in LI 157 are of the same order of magnitude as in the other objects, 
arguing for very similar physical conditions of the molecular gas in all these outflows from 
Class sources. We only note a lower value of the initial ortho-to-para ratio, that is on 
average ~ 1.8, and a significantly higher (3 (up to 5.5). While the former might indicate a 
lower temperature of the gas prior to the passage of the shock, the latter might not reflect 
intrinsic shock conditions, but rather be due to the LTE assumption under which (3 was 
computed in Paper II. 
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Fig. 12. — Diagnostic diagram of the ratio HD R(3) /R(4) against H2 S(2)/S(0), as a function 
of a grid of temperature (200 < T < 400 K, continuos lines) and density (10 4 5 < n(H 2 ) < 10 6 
cm -3 , dotted lines). The effect of introducing atomic hydrogen as a collision partner is shown 
with a long-dashed curve, for \og(n(H 2 ))=4:.5 and /=0.5: this curve coincides with the curve 
with \og[n(H 2 )}=4:.7 and f=0. Hence, an increase of / implies lower n(H 2 ). 
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Fig. 13.— L1448 images of H 2 S(0)/S(l)/S(2) and HD R(3)/R(4) in a 10 x 20 arcsec box 
around L1448-mm. The same emission scale is used for an easy by-eye comparison of the 
relative brightness. 
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6. Energetics 

According to models for non-dissociative shocks (e.g. Kaufman & Neufeld 1996), H2 
emission represents one of the major channels for radiating away the mechanical energy of 
the flow. For example, for pre-shock densities lower than 10 5 cm" 3 and shock velocities 
larger than 20 km s _1 , the fraction of radiation in form of H 2 lines is between 40% and 70%. 
The luminosity of H 2 S(0)-S(7) rotational lines for the three outflows presented here has 
been estimated in Paper I (their Table 3). Other contributions from bright H 2 lines can be 
evalutated by considering detailed predictions of shock models: for example Smith (1995), 
predicts at least ten near- infrared lines (mainly coming from the v=l level) with intensity 
comparable to that of the 2.12 //m; since from our images we measure L(2.12) ~ 0.1-0.2 
L(0-0), we can roughly estimate the total H 2 luminosity by simply doubling that of the 
0-0 lines. For example, we get in this way L(H 2 ) determinations for L1448 and BHR71 in 
agreement within a factor of two with those estimated by Caratti o Garatti (2006) in LTE 
approximation. 

Other important contributions to the total cooling come from species emitting in the 
far-infrared (Kaufman & Neufeld 1996). These have been measured in shocks from Class 
sources by Giannini, Nisini, & Lorenzetti (2001); basing on ISO observations they found 
that the FIR cooling (L(FIR)=L(OI)+L(CO)+L(H 2 0)+L(OH)) is between 10~ 3 and 10" 1 
L , the dominant cooling species being CO and H 2 0. 

In L1448, the L(0-0) luminosity of the northern lobe is 0.029 L , i.e. more than twice in 
the southern lobe (recall, however, that the most distant bow-shock of the southern lobe was 
not mapped). This difference can be partially due to the larger extension of the northern lobe 
(by about a factor 1.5), but likely also reflects the presence of internal shocks caused by the 
interaction of the two outflows coming from L1448-mm and the triple system L1448-IRS3. 
This is also demonstrated by the higher densities measured, on average, in the northern 



Table 3. Line Cooling 1 . 



Object 


L(0-0) 


L(H 2 ) 


LcooZ 


Lfcin 


L coo i/Lf)o! 




(L© ) 


(L©) 


(L©) 


(L©) 




L1448 


0.044 


0.09 


0.44 


0.33 


0.022 


BHR71 


0.048 


0.096 


0.096-0.36 




0.0005-0.02"; 0.06-0. 2 b 


NGC2071 


~1.14 


2.3 


2.3-10 


1.6-180 


0.0004 -0.02 c 



1 References are given in the text 

a : computed for IRS1; ^computed for IRS2; c :value referring to the whole 
cluster of 8 sources. 
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lobe, and, in particular, at the peak N: here we register the highest luminosities, also in 
comparison with peak S in the southern lobe (see Table CP . The higher luminosity of the 
northern lobe is evident also in the FIR cooling, where L(FIR) is 0.076 L (against 0.047 L Q 
in the southern lobe), even if the largest value is measured close to the central source where 
L(FIR)=0.23 L (Nisini et al. 2000). 

As shown in Table|3j summing up all the contributions (H 2 + FIR), we get an estimate of 
the total shock cooling, L coo i ~ 0.44 L Q . Since we are not able to distinguish the single 
contributions of each of the two outflows, we compare this quantity with the sum of the 
bolometric luminosities of the two driving sources: these are 11 and 9 L for L1448-IRS3 and 
L1448-mm, respectively (Andre et al. 2000). The so derived L coo i/L bo i ratio is in agreement 
with other determinations for Class sources (Giannini, Nisini, & Lorenzetti, 2001). L coo i is 
also comparable with the most recent estimate of the kinetic energy of the molecular outflow 
(Curtis et al. 2010). Hence, the energy released in the multiple shocks inside the jet is shared 
out roughly in equal parts between line radiation and mechanical energy of the outflow. 

In the case of BHR71, we have computed separately the L(0-0) luminosity for the 
northern lobes of the two outflows driven by the two sources IRS1 and IRS2. These are 
L(0-0)= 0.019 L and 0.008 L , respectively. The two southern lobes cannot be separated 
observationally; if we assume the same proportions as for the northern lobes we get for the 
IRS1 (IRS2) outflow L(0-0)= 0.015 (0.006) L . With the same prescriptions adopted for 
L1448, the rough estimate of the total H2 cooling for the two outflows is 0.068 L and 0.028 
L . Since there are no estimates for the FIR cooling, we assume L coo ; ranging between L(H 2 ) 
and 4x L(H 2 ), under the assumption that H 2 makes at least a 25% contribution to L coo ;. 
With the new determinations of the bolometric luminosities of the two sources, based on 
interferometric observations, of 13.5 and 0.5 L (Chen et al. 2008), we derive L coo //L;, /=5 
10~ 3 -2 10~ 2 and 5.6 10~ 2 -2.2 10" 1 for IRS1 and IRS2, respectively, (see Tabled . 

Our estimate of the H 2 cooling in NGC2071 is certainly underestimated, since our 
observations cover the south lobe only partially. The separate computation of L(0-0) in the 
two lobes gives L(0-0)=0.57 L and 0.19 L for the northern and southern lobe, respectively. 
If we attribute the same luminosity for each lobe, we get L(0-0)~1.14 L , and L(H 2 )~2.3 
L . Again, if we assume that H 2 contributes at least 25% to the total cooling, a reasonable 
range for L cooi is between 2.3 and 10 L . We compare this rough evaluation with the 520 
L total luminosity of the cluster, composed of eight distinct infrared sources (Butner et 
al. 1990, Walther et al. 1993). The ratio L coo i/Lboi is within the typical range observed 
for outflows coming from Class sources (see Tabled]). Finally, our estimate of L cooi falls 
between the two more recent estimates of the kinetic luminosity of the high velocity gas of 
1.6 L (Stojimirovic, Snell, & Narayanan 2008) and 180 L (Buckle et al. 2010). These 
estimates are so discrepant that it is impossible to perform a meaningful comparison with 
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our L coo i determination. 

7. Conclusions 

We have presented the H 2 pure rotational lines from S(0) to S(7) observed with IRS- 
Spitzer in three outflows from Class protostars: L1448, BHR71, and NGC2071. We have 
analyzed the data by assuming a gas temperature stratification that follows a power law 
distribution. The mid-IR data have also been combined with near-infrared H 2 rovibrational 
lines and mid-IR HD lines detected in few spots to derive the physical parameters under a 
NLTE approach. Our main results are summarized in the following: 

- We have obtained detailed maps of gas column density, OPR and temperature assuming 
the LTE approximation. All the flows present similar characteristics in regard to the 
covered ranges of these parameters, and the derived values are similar to the same 
quantities measured in Paper II in LI 157. Typically the column density ranges between 
1-7 10 20 cm' 2 , corresponding to a mass of gas of the order of 1CT 2 M . The only 
exception is NGC2071, which is definitely more massive than the other two (M~ 0.3 
M ). 

- A large range of temperatures, between few hundreds to few thousands of Kelvin, 
contribute to the H 2 emission, in all three of the flows. The hottest components are 
found typically along the jet axis and in warm spots, corresponding to Herbig-Haro 
objects or to small regions where HD emission is also detected. 

- Significant deviations from the equilibrium value of the ortho-to-para ratio are evident 
in all the three outflows. On average we find values around 2.2, with a large range that 
goes from less than one to three. However, some trends can be identified, with the 
northern lobes of L1448 and BHR71 closer to the equilibrium than in the south, and 
with a decrease in NGC2071 going from the internal toward the apexes of the flow. 
Generally speaking, higher OPRs are found at the location of temperature peaks, even 
if low OPR values can sometimes be found also in regions with relatively high temper- 
atures, as for example in the southern lobe of BHR71. These apparently contrasting 
results might be explained if the time elapsed from the last shock event were no long 
enough to allow reaching the low OPR values. 

- H 2 rovibrational lines, and in particular the 1-0 S(l) line at 2.12 jum, have been used 
in conjunction to the rotational lines to relax the LTE approximation, and to derive 
detailed maps of the post-shock molecular gas density and the power law index that 
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regulates the temperature stratification. As far as this latter parameter is concerned, 
we find values always larger than 3.8, in agreement with shock models for multiple 
C-types bow shocks with apex velocities able to partially dissociate H 2 . 

- H 2 densities are roughly between 10 6 -10 7 cm -3 ; these values, however, can decrease by 
up to a factor of ten if a significant fraction of atomic hydrogen exists; moreover, from 
the analysis of the HD emission we demonstrate that indeed a density stratification does 
exists, with the lowest densities corresponding to the lowest temperature components. 

- In a few selected regions we were able to construct rotational diagrams up to excitation 
energies of ~ 20 000 K and v < 3. From these diagrams we are able to estimate the 
time needed for the para-to-ortho conversion as a function of the fraction of the atomic 
hydrogen. In HH320 and HH321, in particular, we used the results of recent shock 
models to estimate that fraction, finding that hydrogen is almost all in molecular form. 
This result, although regarding small regions inside the flow, gives a different picture 
of that found in LI 157, where the fraction of atomic hydrogen has been found to be ~ 
0.1-0.3. 

- We have studied in detail a region, about 20" in length, close to the L1448-mm source. 
Here very different conditions seem to occur in the two lobes: in the northern one all 
the H2 rotational lines are detected, indicating a temperature ~ 300 K; conversely in 
the south lobe only the fundamental line is detected, along with the HD lines R(3) 
and R(4). These findings are consistent with a gas temperature lower than 200 K and 
a local enhancement of the column density: we expect consequently that a different 
chemistry could have taken place in the two lobes. 

- An estimate of the radiative cooling of H2 was derived and compared, when available, 
with the cooling of abundant molecules emitting in the far-infrared. The estimated 
total cooling, if compared with the kinetic energy of the swept-out outflow, supports 
the scenario in which the shock from which H 2 is emitted is also capable of accelerating 
the molecular outflow. 

Finally, we note that the objects analyzed in the present paper are going to be observed 
with the far-infrared and submillimetre telescope Herschel: the map parameters we have 
provided here will be essential for the interpretation of the forthcoming data, and will be 
helpful for deriving a more unified view of the physical and chemical processes governing the 
ongoing shocks. 
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